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Abstract 

The effects of binary tidal forces on transport within an accretion disk are studied 
with a time-dependent hydrodynamical model of a two-dimensional isothermal ac- 
cretion disk. Tidal forces quickly truncate the accretion disk to radii of order half 
the average radius of the Roche lobe, and excite a two-armed spiral wave that re- 
mains stationary in the rotating reference frame of the binary system. We measure 
an effective a of order 0.1 near the outer edge of the disk in all of our models, 
independent of the mass ratio, Mach number, and radial density profile. However, 
in cold disks with high Mach number, the effective a drops rapidly with decreasing 
radius such that it falls below our threshold of measurement (~ 10 -3 ) at a radius 
of only one third the tidal radius. In warmer disks where the Mach numbers remain 
below 20, we can measure an effective a down to radii 10 times smaller than the 
maximum size of the disk. 
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1 Introduction 



Accretion disks play a prominent role in modern astrophysics, yet we remain 
largely ignorant of their detailed structure and local dynamics. The impor- 
tance of angular momentum transport has been understood for a long time 
(Shakura & Sunyaev 1973, Lynden-Bell & Pringle 1974), but for many years 
our understanding of disks has been built on a parameterization of some un- 
known transport mechanism. In recent years much progress has been made in 
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elucidating the role of magnetohydrodynamical turbulence in mediating an- 
gular momentum transport in accretion disks. A recent review of our current 
understanding of transport is given by Balbus & Hawley (1998). 

In addition to local sources of transport such as MHD turbulence, global pro- 
cesses may also contribute to the transport of angular momentum in accretion 
disks. In particular, spiral shock waves can act as a local sink for angular mo- 
mentum (Spruit 1987, Larson 1989). Because these waves are traveling slower 
than the local Keplerian velocity, they contain negative angular momentum, 
and any dissipation at the shock will decrease the angular momentum of the 
orbiting gas. While spiral waves can be produced by a variety of sources, 
they arise naturally in binary star systems such as X-ray binaries, cataclysmic 
variables, and binary proto-stars, where the gravitational pull of a binary com- 
panion creates a two-armed spiral shock wave that co-rotates with the binary 
system. 

Such spiral waves have been observed in several numerical studies, beginning 
with the early work of Sawada et al. (1986) and continuing with recent at- 
tempts to model the observed spiral pattern in IP Peg (Godon et al. 1998, Ar- 
mitage & Murray 1998). Despite this large body of numerical work, there has 
been relatively little progress in the quantitative measurement of the effective 
transport of spiral shocks in hydrodynamic disks. This is due in large part to 
the fact that it is easy to set up a hydrodynamic simulation of a Keplerian 
accretion disk, but it is very hard to do it accurately enough to measure the 
relatively small transport expected due to spiral waves. One needs high nu- 
merical accuracy, a large range in radii, and long simulations covering several 
orbital periods. Rozyczka & Spruit (1993) ran a simulation without a tidal 
stream for tens of orbits. While their models did exhibit a stationary two- 
armed spiral pattern, the shock dissipation was not enough to keep up with 
the cooling in their model. As the gas in the disk cooled and the Mach num- 
ber increased, the spiral shocks became weaker and eventually faded away. 
Savonije et al. (1994) performed careful numerical simulations to study tidal 
effects, making a quantitative comparison with analytic calculations of the 
linear response of the disk to tidal forces. They found that spiral shocks were 
much diminished in cold disks with Mach numbers of order 25 or higher, a 
feature they attributed to a mismatch between the wavelength of the tidal re- 
sponse of the disk and the length-scale of the driving tidal force. Godon (1997) 
ran relatively high resolution simulations using a markedly different numeri- 
cal technique (a hybrid spectral method) and found results similar to previous 
authors, namely a strong two-armed spiral wave. Note that both Savonije et 
al. (1994) and Godon (1997) had disks that covered only a limited range in 
radii (4 and 5 respectively), and so could not follow the propagation of spiral 
waves deep into the disk. 

The results of these simulations suggest that spiral waves are a robust feature 
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of accretion disks in binary systems, and that these spiral shocks can indeed 
transport mass and angular momentum through the disk. They do not, how- 
ever, provide a quantitative measure of this transport. The goal of this paper 
is to measure an effective a due to spiral shocks excited in accretion disks by 
the tidal forces of a binary companion. 



2 An Idealized Accretion Disk 

We begin by defining an idealized theoretical model of an accretion disk that 
we can simulate with a standard time-dependent hydrodynamics code. Of 
particular concern is the treatment of the thermal structure of a disk, which is 
determined by a balance between local energy deposition and vertical energy 
transport via convection and radiative diffusion. If an adiabatic equation of 
state is used and no outlet is provided to remove the thermal energy, the 
disk will gradually heat up. This approach has been used by several authors 
(Sawada et al. 1986), but the increasing sound speed changes the model as the 
simulation evolves, and ultimately leads to hot disks with a correspondingly 
large scale height and small Mach number. 

In order to include the thermal structure of the disk, a model must include 
energy transport in the vertical direction. One way to do this without includ- 
ing radiative diffusion and three dimensions (a difficult prospect even for the 
fastest computers) is to model energy loss by including an ad hoc cooling 
function designed to remove energy at the same rate one would expect due 
to thermal radiation from the disk surface. Rozyczka & Spruit (1993) used 
this approach, specifying a local cooling rate proportional to T 4 . While this 
approach will not necessarily produce a realistic temperature structure in the 
presence of shocks, it does allow the use of a realistic value of 7 while keeping 
the temperature in the disk well below the local virial temperature. Godon 
(1997) avoided the issue of shock heating by assuming a polytropic equation 
of state, P = KYP , where K is a predefined constant. Savonije et al. (1994) 
took an intermediate approach and included an energy equation but added 
local cooling to maintain a polytropic equation of state. 

We take a slightly different approach and assume an isothermal equation of 
state. This is equivalent to the polytropic equation of state in the limit that 7 = 
1. The assumption of an isothermal gas provides for a simpler problem, fewer 
parameterizations, and allows for faster computation. However, the analysis 
of Spruit (1987) suggests that spiral waves will not propagate for a strictly 
isothermal equation of state, so we may be underestimating the role of spiral 
shocks in this work. 

An ideal, isothermal gas evolves according to the Euler equations, which can 
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be written in conservative form as: 

dp 



m + V • (pu) = (1) 
^ + V • (puu) + VP = pa, (2) 



where p is the mass density, u is the velocity, and the pressure is given by 
P = c 2 s p with c s being the isothermal speed of sound. In a coordinate frame 
Co-rotating with the binary system, the acceleration of the gas can be written 

as 

a = -V$ -2fixM 



where the last term is the Coriolis force due to the Co-rotating coordinate 
system, Q is the orbital frequency of the binary system, and the effective 
binary potential is given by 

GMi GM 2 1-. 2 

$ = 77 (n x r cm y. 

r\ r 2 2 



The last term is the centrifugal force due to the Co-rotating system, and r x is 
the distance to the accreting star, r 2 is the distance to the donor star, and r cm 
is the distance to the center of mass of the binary system. We have normalized 
our model such that the unit of distance is given by the binary separation and 
GM\ — 1. In these normalized units the orbital period is P = 2n/y/l + q, 
where q = M 2 /M 1 is the mass ratio. 

We further restrict our idealized model to two dimensions by assuming no 
vertical motions nor variations with height in the disk. When using a cylin- 
drical coordinate system, this is equivalent to replacing the density in the 
disk with a surface density, S = 2pH, and assuming a constant scale height, 
#(Godon 1997). 

Because our objective is to measure the tidal effects on transport in an accre- 
tion disk in a binary system, we need some quantitative measure of transport 
through the disk. The most direct physical quantity we could measure is the 
mass accretion rate, given by M = 2nrYiV r in two dimensions. The surface 
density has no particular scale (we chose E = 1 in our flat disk models); the 
important parameter is really the radial velocity, v r . However, rather than 
reporting the average radial velocity (as compared to, say, the sound speed), 
we have followed the historical convention of this subject and reported the 
effective value of a. Shakura & Sunyaev (1973) introduced a as a means of 
parameterizing some unknown form of angular momentum transport, presum- 
ably related to local turbulence in some fashion. While the transport being 
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measured in our models is a global phenomenon rather than a local process, 
we can still refer to an effective a that would give the same mass accretion 
rate in an a-disk model that we find in our simulations. 

The standard a-disk theory gives the mass accretion rate in terms of a as 
M = 37rac s HH, where H is the local scale height in the disk. Using the fact 
that H = Rcs/vtj, for a cold, thin disk, we find an equation for a: 

2v lP (v r ) 

where (v r ) is a density-weighted average of the radial velocity: 

_ J 27r Y,v r d<f) 

{Vr) ~ • 

Thus, to measure a we need to keep track of the local radial mass flux (f r S) 
and the local surface density. 

It is important to point out that our ideal accretion disk is not expected to be 
in strict steady state. We have no a priori reason to expect the mass accretion 
rate to be independent of radius. Mass accretion in this model is driven solely 
by the action of spiral shock waves, whose strength will certainly depend on 
the radius within the disk. The propagation of these shocks will also depend 
on the density profile in the disk, but there is no guarantee that a combination 
of density profile and shock propagation exists that leads to a constant mass 
accretion rate. Furthermore, even if there did exist a steady-state solution, we 
may not, in practice, be able to evolve a simulation long enough to reach a 
true steady state. 



3 Computational Method 

Our idealized accretion disk model is evolved using a modified, isothermal ver- 
sion of the hydrodynamics code VH-lhttp : //wonka. physics .ncsu. edu/pub/VH-1, 
run in two-dimensional cylindrical coordinates (r, 0) centered on the accret- 
ing star. This code is based on the Lagrange-remap version of the piece-wise 
parabolic method (Colella & Woodward 1984). The isothermal version is sig- 
nificantly faster than the adiabatic code because there is no need to solve an 
energy equation and there exists an analytic solution to the Riemann problem, 
avoiding a costly iterative solution as required in an adiabatic code. 

There are essentially two modifications to the standard distribution of the 
isothermal version of VH-1: Forces to account for the rotating binary potential 
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and changes in the remap step to allow conservation of angular momentum. 
This second step is vital to any quantitative study of accretion disks, for it is 
the transport of angular momentum that allows accretion to take place. If a 
numerical method does not conserve angular momentum, accretion may result 
solely from, or be completely wiped out by, numerical error. 

Conserving angular momentum in VH-1 is a fairly simple task. When comput- 
ing the hydrodynamic evolution in the direction the radius is held constant 
so conserving angular momentum is equivalent to conserving momentum, 
which VH-1 already does. When computing the evolution in the radial direc- 
tion the angular momentum in a zone is left unchanged during the Lagrangian 
hydrodynamic evolution. This leaves the remap step in the radial direction as 
the only place in the standard code where angular momentum is not conserved. 
Fortunately, conserving angular momentum in this step is trivial: instead of 
interpolating on and remapping momentum, one need only interpolate on 
rVcj) and remap angular momentum. Doing so will conserve the total angular 
momentum on the grid to within machine roundoff error (if boundaries are 
reflecting and forces are only radial). 

An additional concern is the local diffusion of angular momentum. Even if 
angular momentum is conserved globally, it may diffuse across the numerical 
grid, producing anomalous transport. Savonije et al. (1994) modified their code 
to interpolate on the quantity Xr 1 / 2 -^ (which is flat in a Keplerian disk with 
constant E) in order to minimize the diffusion of angular momentum. This 
prevented unwanted instabilities from developing in their disk, although it did 
not necessarily guarantee conservation of angular momentum. Fortunately, 
VH-1 exhibits very little diffusion of angular momentum in this problem, a 
feature we attribute to the higher-order interpolations used in the piece-wise 
parabolic method. 

The Lagrange-remap method used in VH-1 allows a simple, accurate mea- 
sure of the radial mass flux need to compute an effective a. Advection in 
VH-1 is treated in the remap step, where gas that has been evolved on a 
Lagrangian grid is remapped back onto a fixed Eulerian grid. For example, 
if the Lagrangian interface between zones 1 and 2 moves to the right during 
a timestep, the remap step takes the mass (and momentum) within a sub- 
shell defined by the original zone interface and the final Lagrangian interface 
and instantaneously moves it from zone 1 (where this mass was before the 
timestep) to zone 2 (where it has advected during the timestep). The mass 
flux through the Eulerian interface during the timestep is just the mass in the 
sub-shell that was remapped from zone 1 to zone 2. So, to compute an average 
mass flux, t> r S, we summed up all the mass being remapped across a given 
radial boundary (i.e., for all values of 0) during a prescribed interval of time 
(typically 50 time steps). Dividing by the average (over time and 0) density 
in these sub-shells, we can calculate an effective value of a. 
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The initial conditions for the disk simulations include a low-density hole in 
the first 5 radial zones (10~ 3 of the disk density) to act as a buffer between 
the inner edge of the disk and the inner boundary, zero radial velocity, and 
= r -1 / 2 — £lr. The last term in accounts for the centrifugal force in 
the rotating reference frame, and provides for an initial disk much closer to 
equilibrium than pure Keplerian rotation. 

The boundaries are periodic and the radial boundaries are constructed to 
allow material to flow off the grid, but not onto the grid. At the inner radial 
boundary the velocity is set to zero in the ghost zones (extra zones used 
for interpolation on the real grid), creating an imbalance of forces at the 
inner edge of the hydrodynamic grid that drives quick accretion off the inner 
edge. Contrary to some previous numerical work (Matsuda et al. 1990), the 
exact boundary conditions at the inner edge do not significantly alter the 
simulations. Even with a reflecting boundary, material piles up at the inner 
edge of the grid without influencing the spiral shocks at larger radii. The outer 
boundary also does not exert much influence on the long term evolution of the 
disk. Once tidal forces trim down the outer edge of the disk, the density in the 
outer zones is very small so there is negligible flux of mass and momentum off 
the outer edge of the grid. 

Most of our simulations were run on grids with a radial span of a factor of 20, 
with the radial zone size gradually increasing with increasing radius to provide 
consistent resolution over the full range of radii in the simulation. The outer 
radius was chosen to be approximately the distance to the inner Lagrange 
point, allowing plenty of room between the outer edge of the tidally-truncated 
disk and the outer grid boundary. One would like to make the inner edge of 
the numerical grid as small as possible, but restrictions on the timestep in an 
explicit hydrodynamic code make small inner radii very expensive to compute. 
For stability, the timestep must be less than rA0/t^. But = r -1 / 2 in an 
accretion disk, so the limiting timestep decreases as r 3 / 2 . Evolving a disk that 
extends down a factor of 20 in radius thus requires a small timestep and a 
corresponding large number of cycles; some of our simulations ran for over 10 6 
cycles. The choice of 20 for the radial span was based on the radial decay of 
a in our standard model, as we will see in the following section. This inner 
radius also corresponds roughly to the point of closest approach for a ballistic 
stream of gas originating at the L\ point in the binary system. Since our goal 
in future work is to quantify the effects of the tidal stream on transport within 
the disk, we have extended our disks in this work down to this small radius. 

The issue of numerical resolution was approached empirically. Keeping in mind 
that our goal was a quantitative measure of transport effects caused by tidally- 
induced spiral waves, we ran several simulations with different numerical res- 
olutions. We found relatively little difference in the density structure of the 
disk and the computed value of a when we used 115 zones in the direction 
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Fig. 1. Resolution of spiral waves as a function of the number of radial zones, n r . 
The grey scale represents gas density, which reaches a maximum of approximately 
3 times the initial disk density in the spiral waves. 

compared with 200 zones, so we opted for the former in most simulations. For 
the radial direction we found significant differences between runs with 115 and 
200 radial zones, but relatively little difference between runs with 200 and 300 
radial zones. Figure 1 shows the loss of definition in the spiral waves when 
the radial resolution is too low. The computed value of a in the run with 200 
zones evolved very similar to that in the run with 300 zones, but the run with 
100 zones looked completely different (and in fact was negative throughout 
the evolution). More radial zones were required for the colder disk simulations 
described in Section 4.2 

Finally, a small amount of dissipation was included in the code in the form of 
wiggling the grid back and forth a small amount in the radial direction (Colella 
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& Woodward 1984). This was necessary in order to suppress two different 
numerical instabilities found in these disk simulations. The first creates a 
zone to zone striping in the radial direction, and the second leads to epicyclic 
waves, usually (but not always) appearing at the inner boundary of the disk. 
The first problem is related to the conservation of angular momentum. In 
a one-dimensional (radial) problem using Lagrangian coordinates, if a given 
zone has too much angular momentum it will drift outwards. However, this 
motion will drag out neighboring zones, which will find themselves having too 
little angular momentum for their new radius and will drift inwards. This can 
lead to a striping effect with the radial velocity in each zone changing sign 
from zone to zone. We were not able to fully understand the second problem, 
but we found that it was worse (in the absence of grid wiggling) if the ratio 
of the sound crossing time across a radial zone to the orbital period at that 
radius was larger than some critical value. For all practical purposes these 
instabilities were eliminated by the grid wiggling, which in turn had negligible 
effect on the structure of the disk or the computed value of a. 

As a test of our computational method we ran a simulation without the binary 
potential, i.e., with only the central force of the accreting star's gravity. We 
instituted reflecting boundary conditions at both the inner and outer radii, 
so the total mass and angular momentum on the grid should remain fixed. 
To make the test meaningful, we perturbed in a small band of radii in the 
middle of the grid to create wave motions. The mass and angular momentum 
did remain constant to within machine round off error. The local radial veloc- 
ities remained less than ~ 10% of the sound speed, and the computed average 
of a was ~ 1CT 3 . This then represents the limit to which we can measure an 
effective a in our binary simulations. 



4 Numerical Simulations 

This section describes seven hydrodynamic simulations varying the sound 
speed, mass ratio, and density profile in the disk. The parameters describ- 
ing each of these models are listed in Table 1. 

4-1 Standard Model 

For the purposes of presenting our results we have adopted a standard model 
in which q = 1 and c s = 0.25. We will describe this model in detail, and then 
go on to discuss how other models differ from this standard. With the adopted 
sound speed, the Mach number of the disk material at the inner edge of the 
disk is ~ 28, comparable to some of the best previous numerical work, but still 
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Table 1 

Input parameters (sound speed c s , mass ratio q, and surface density index n) for 
the hydrodynamic simulations. The tidal radius, Rt, and the Mach number of the 
Keplerian flow at half the tidal radius are given, as well as the ratio of the tidal 
radius to the radius of the Roche lobe of the accreting star. 



Model 


c s 


q 


n 


R T 


Rt/Rl 


M(R T /2) 


standard 


0.25 


1.0 





.214 


0.56 


12.2 


colder 


0.125 


1.0 





.264 


0.69 


22.0 


coldest 


0.0625 


1.0 





.300 


0.79 


41.3 


steeper 


0.25 


1.0 


1 


* 


* 


12.2 


steepest 


0.25 


1.0 


3 


* 


* 


12.2 


lowmass 


0.25 


0.2 





.307 


0.59 


10.2 


highmass 


0.25 


5.0 





.185 


0.74 


13.2 



lower than one might expect in real disks. Tidal forces quickly truncate the 
disk at an outer radius of 0.21, which is roughly half of the averaged Roche 
lobe radius as given by the analytic expression of Eggleton (1983). This disk is 
significantly smaller than predicted by theoretical arguments for infinitely cold 
disks(Paczynski 1977), but one expects that finite pressure effects will increase 
the effectiveness of tidal forces in truncating the disk. At this tidally-truncated 
outer edge the disk material orbits with a Mach number of ~ 8. 

The images and associated animations in Figure 2 show the strong spiral struc- 
ture that develops in our standard model. The two-armed spiral wave remains 
relatively fixed in the co-rotating reference frame of the hydrodynamic simu- 
lation. The shape of the spiral is determined by the inward radial propagation 
of the wave at a velocity of order the sound speed and the shearing of this 
wave by the orbital motion. The result is a trailing angle (between the spi- 
ral shock and a circle of given radius) given by tan0 m (c s /Qr) (Godon et 
al. 1998), or the inverse of the Mach number of the orbital motion, M. Larson 
(1989) suggested a radial wave speed of ~ 1.5c s , whereas our simulations show 
a spiral shape corresponding to a wave speed of ~ 1.3c s . The density jump at 
the spiral shocks varies from a maximum of ~ 2.1 near the outer edge of the 
disk down to a value of ~ 1.4 near the inner edge. The corresponding shock 
Mach numbers of order ~ 1.3 (given by the square root of the compression) 
are in agreement with the low Mach numbers predicted by Spruit (1987) and 
with the predicted shape of the spiral pattern. Presumably the strength of the 
spiral shocks would continue to decrease if we had extended our simulation 
down to even smaller radii, as the spiral shock becomes wound tighter and 
tighter with increasing orbital velocity. 

The initial response of the disk to the binary potential and the subsequent 
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Fig. 2. Density in the standard model at time t = 9. The left image shows the entire 
grid spanning the Roche lobe of the accreting star, while the right image shows 
the central region magnified by a factor of 5. The binary companion is to the right 
of the image, and the disk orbits in a clockwise sense. Black corresponds to high 
density gas, with saturation occurring at a density of 3. The black circle in the right 
image marks the inner edge of the numerical grid at r = 0.02. The accompanying 
animations show the evolution of the disk from the same two perspectives. 




time 



Fig. 3. Time evolution of a in the standard model at a radius of r 
Keplerian period at this radius is Pk ep ~ 0.105. 



0.065. The 



quasi-steady transport are illustrated in Figure 3, which plots the effective 
a at a radius of r 0.065. There are very large waves generated by our 
crude initial conditions, as shown by the large variations in a shortly after 
t — 0, but these quickly dissappear on a timescale of less than one rotation 
of the outer disk edge. The behavior of a does not appear to settle into a 
quasi-steady behavior until after approximately one orbital period (4.44 in 
this model), but from 1 to 3 orbital periods the behavior appears unchanged. 
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Fig. 4. Time-averaged radial profile of a in the standard model for times of 7-9 
(solid), 9-11 (long dash), and 11-13 (short dash). The dotted line shows a at late 
times (t = 25 — 27) after other wave modes have appeared in the outer regions of 
the disk. 

The high frequency oscillations clearly visible in Figure 3 correspond to the 
local Keplerian frequency (Pk ep = 2nr 3 ^ 2 w 0.105). 

The variation of a with radius is shown in Figure 4. The value of a is as high 
as 0.1 near the outer edge of the disk where strong spiral shocks are driven by 
tidal forces, but appears to gradually decline with time at large radii. Inside 
of r ~ 0.1, a appears relatively constant with time but drops with decreasing 
radius. Although a(r) is not well fit by a power-law, we can approximate a in 
the central region of the disk with a r 1A . At late times other wave modes 
appear (see the animation in Fig. 2), slightly reducing a at large radii and 
enhancing it at small radii. The result is a more uniform value of a throughout 
the disk. However, more analysis of these other wave modes is needed before 
one can be confident of their role in angular momentum transport. For this 
reason we concentrate our attention on the two-armed spiral waves. 

These values of a due to the two-armed spiral shocks are consistent with the 
shock strengths quoted above, but the radial dependence is much steeper than 
the predictions of Larson (1989), and the absolute values are much higher than 
the analytic predictions. For a spiral shock with Mach number M s , Larson 
(1989) derives an effective a in an isothermal disk of 



For an isothermal shock, M 2 is equal to the compression ratio. The range in 
M 2 S from 1.4 to 2.15 gives a range in a of 0.005 to 0.1, consistent with the 



a w 0.07(M S 2 - l) 3 . 



(4) 
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Fig. 5. Radial profile of £ (averaged over (f>) in the standard model at various times 
beginning with t = 9 (solid line) and ending with t = 16 (dotted line). 

data plotted in Figure 4. 

Larson (1989) goes on to predict that the strength of tightly wound spiral 
shocks in a cold, thin disk should be proportional to (cos#) 3 / 2 , where 9 is the 
trailing angle of the spiral wave. Estimating cos 9 ~ 1.5c s /v , Larson (1989) 
finds an effective a of 



The prediction for an isothermal disk with c s = 0.25 is a ~ 3.2 x 10 3 r 3 / 4 , 
orders of magnitude below that found in our simulations. 

A simple estimate of the timescale for gas to move through the disk is r/v r , 
which for this model ranges from 60 at the outer edge of the disk to ~ 300 
near the inner edge. Having evolved our simulation for a time of only 13, we do 
not expect our model to have settled into an equilibrium. We can in fact see 
small changes in the density profile of the disk as plotted in Figure 5. Despite 
the fact that the density is slowly increasing over most of the disk, the total 
mass on the computational grid has decreased by almost 1% during the time 
from t — 9 to t — 16. 

The two-armed spiral pattern seen in Figure 2 is not the only wave present in 
the disk. On a timescale of several binary orbits other wave modes begin to 
appear. By viewing the animation accompanying Figure 2, one can identify a 
purely radial wave (m = 0) and a precessing eccentric wave (m = 1) growing 
in the vicinity of the outer disk edge. While these waves continue to grow to 
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some apparent saturation level, they do not completely disrupt the two-armed 
spiral shock, and a remains at a significant level 4. Similar wave modes have 
been seen in the simulations by Godon (1997). We will discuss them further 
in the next section, for they are more apparent in colder disks. 



4-2 Dependence on Mach Number 

The sound speed in our standard model was chosen to produce a significant 
value of a, not because it was astrophysically interesting. In most accretion- 
powered binary systems the Mach number of the orbital motion in the disk 
will be much higher than we have assumed here. A colder disk is both more 
challenging to simulate and more difficult to measure the expected smaller a. 
The effect of spiral shocks on transport through an accretion disk diminishes 
rapidly with increasing Mach number because the spiral shocks are wound 
tighter, making the relative shock Mach number smaller, leading to less dissi- 
pation and effective transport. 

In addition to the challenge of measuring small effects, we were faced with 
several numerical problems when simulating colder disks. The first problem is 
one of resolution. Because the spiral arms are wound tighter in colder disks, 
we found that to accurately resolve the oblique shocks of the spiral arms, we 
needed higher resolution in the radial direction. Another difficulty plaguing our 
cold disk simulations was the increase in numerical " noise" with higher Mach 
number. In one-dimensional (radial) simulations of Keplerian disks excited 
with radial waves, we found a monotonic increase in the average radial velocity 
with increasing Mach number in the disk. While the effective a from this 
numerical noise remained low (a = 4x 1CT 3 for c s = .0125), it could affect the 
measurement of transport in very cold disks. Finally, the numerical instabilities 
described in the appendix are more apparent at higher Mach numbers. 

To investigate the influence of sound speed on our idealized accretion disks we 
ran simulations similar to our standard model but with sound speeds of 0.125 
(with 300 radial zones) and 0.0625 (with 400 radial zones). As expected, the 
tidal forces have less of an effect on the colder disks. As seen in Figure 6, the 
spiral density waves are wound tighter and fade faster with decreasing radii 
in the colder, higher Mach number disks. The two-armed spiral has almost 
completely disappeared by the inner edge of the disk with c s = 0.125, and is 
present only beyond r 0.08 in the coldest disk. 

A one-armed spiral is clearly visible in the interior of the coldest disk. This 
wave originates in the outer regions of the disk and propagates inward, decay- 
ing away below a radius of r pa .04. There is no trace of this wave at the inner 
edge of the disk at any time during the simulation. A similar wave can be seen 
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in most of the other simulations, showing up either as an initial transient, 
or as a growing mode at late times in the simulation. However, this wave is 
much harder to identify in the other models, both because the single spiral is 
weaker, and also because the double spiral is much stronger in warmer disks 
and hence obscures the single spiral. The one-armed spiral is easily seen in the 
animation of the cold disk with c s = .125, but it is a transient phenomena, 
presumably excited by the initial conditions. At later times there is a faint 
hint of the one-armed spiral, but it is much weaker than the double spiral 
wave. 

The source of this single spiral wave (ignoring for the moment any transient 
created by the initial conditions) may be related to the fact that the double 
spiral wave at the outer edge of the disk is not steady, and in fact is much 
less steady in the coldest disk where the single spiral is most prevalent. In 
contrast to the standard model, there appears to be at least two competing 
waves excited at the tidal radius in the coldest disk model, for the two-armed 
spiral is never steady. 

Godon (1997) also observed a single spiral in his Model 1, which has Mach 
numbers comparable to our coldest disk. He attributed this wave to "viscous 
oscillations" excited by the tidal potential. In contrast to our simulations how- 
ever, Godon (1997) found that the single spiral wave was trapped between the 
reflecting inner boundary and the outer edge of the disk. As a consequence, 
this wave grew unbounded in his simulation and eventually dominated the 
entire disk. In addition, the presence of this wave was closely tied to his choice 
of a Keplerian, reflecting inner boundary. Any deviation from this assumption 
and the single spiral did not appear. In our simulations, which covered a ra- 
dial range four times larger, the single spiral decayed away before reaching the 
inner edge of the disk, and its presence was unaffected by the choice of inner 
boundary conditions. 

These simulations also provide evidence that tidal truncation becomes less 
effective with increasing Mach number. The tidal radius of the coldest disk 
is ~ 80% of the averaged Roche lobe radius, more in line with the analytic 
predictions of Paczynski (1977). By searching for the largest non-intersecting 
periodic orbits, Paczynski (1977) found an upper limit (in the absence of 
pressure effects) to the size of an accretion disk. For a mass ratio of unity, 
he predicted a maximum disk radius of 0.3, in good agreement with the tidal 
radius in our coldest model. The warmest model had a tidal radius of 0.214, 
some 30% smaller, leaving the disk radius only 60% of the Roche lobe radius. 

Figure 7 shows the precipitous drop in a as a function of radius for cold disks. 
Analytic arguments suggest that a should drop as c 3 / 2 , but maintain the same 
r 3//4 radial dependence in an isothermal disk. We find that the maximum value 
of a, found near the outer edge of the disk, remains roughly independent of 
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Fig. 6. Snapshots of the disk surface density in models with a sound speed of 
c s = 0.25 (left), 0.125 (center), 0.0625 (right). The accompanying animations il- 
lustrate the increasingly unsteady behavior as the sound speed is reduced. 




Fig. 7. Dependence of the effective a on the sound speed in the disk. The solid 
line corresponds to the standard model with c s = 0.25, the long-dashed line to the 
model with c s = 0.125, and the short-dashed line to the model with c s = 0.0625. 

sound speed, while the radial dependence of a(r) steepens with decreasing 
sound speed. In fact, because the radial decay of a is so steep, we found 
steady mass accretion only for radii above r 0.1 in the coldest disk with an 
outer Mach number of ~ 32. 



4-3 Dependence on Density Profile 



The results of our models so far are characterized by a relatively steady mass 
accretion rate that varies with radius in the disk. This implies that the density 
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Fig. 8. Dependence of the effective a on the density profile in the disk. The solid 
line corresponds to the model with n = 3, the short-dashed line to the model with 
n = 1, and the short-dashed line to the standard model with n = 0. 

profile is changing with time as mass piles up in the disk. Left to evolve for 
a much longer time, we would expect the density profile to steepen, with 
higher densities at smaller radii where the spiral waves are less effective at 
transporting angular momentum. Would this different density profile affect 
the computed value of a in the model? Will the spiral waves weaken as they 
propagate down into regions of higher density? 

We ran additional simulations for several density profiles E oc r _n , and show 
results for n — 0, 1, and 3 in Figure 8. For comparison, an a-disk model for an 
isothermal equation of state yields n = 1.5, and for an adiabatic gas n = 0.75 
(Shakura & Sunyaev 1973). Note however that the gas density scales with 
a steeper function of radius than the surface density. In a three-dimensional 
isothermal model the scale height increases as H oc r 3 / 2 , so the density in the 
equatorial plane varies as p oc r~ 3 . 

In all but the most extreme model, the effective transport within the disk is 
very similar. The value of a is down by at most a factor of ~ 2 in the inner 
disk for the model with n = 1 as compared to the standard model, but the 
radial dependence, maximum values, and tidal truncation all appear similar in 
the models with < n < 1.5. In contrast, the spiral waves decay much faster 
in the steepest disk with n = 3. The spiral waves still extend all the way to 
the inner edge of the disk in this model, but the radial velocity perturbations 
and density compressions associated with these waves are lower than found in 
the other models with smaller n. As a consequence of this faster decay, the 
mass accretion rate is still an increasing function with radius, despite the fact 
that the density varies as £ oc r~ 3 . This result suggests that a true steady 
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Fig. 9. Effect of mass ratio, q, on the spiral pattern, with q = 0.2 on the left, q = 1 
in the center, and q = 5 on the right. Each image has been scaled to the size of 
the Roche lobe. The q = 0.2 model was evolved on a grid that went down to much 
smaller radii. The density scale in each image has been normalized to the average 
density in the middle of the disk. 

state may not exist for a two-dimensional isothermal accretion disk in which 
transport is mediated solely by tidally-induced spiral shocks. 



4-4 Dependence on Mass Ratio 

Our last parameterization is the mass ratio, q. In addition to the equal mass 
system studied in the previous subsections, we have run simulations with a 
mass ratio smaller than unity (accretion from a lower mass companion) and 
larger than unity (accretion from a higher mass companion). The spiral pat- 
terns shown in Figure 9 are remarkably similar for different values of q. The 
spiral waves are slightly more open for small q. Since the disk extends out to 
larger radii where the Keplerian velocity is smaller, the Mach number is lower 
and the trailing angle of the spiral is larger. The strength of spiral shocks, as 
measured by the shock compression, is comparable in all three simulations. 

The relative independence of our model on the mass ratio (at least within the 
limited range that we have simulated) is also seen in Figure 10, which shows 
an effective a that does not vary much with q. The only significant difference 
is the location of the outer edge of the disk. But when the tidal truncation 
radius is compared with the radius of the Roche lobe (Eggleton 1983), Rl, 
even the size of the disk does not change much. Accretion from a low-mass 
companion produces a larger disk, but the Roche lobe is larger by almost the 
same amount, leading to a ratio of Rt/Rl comparable to a disk in a binary 
with equal mass stars. Accretion from a high-mass companion produces a 
smaller disk, but one that occupies a larger fraction of the Roche lobe (Table 

1)- 
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Fig. 10. Dependence of the effective a on the mass ratio, q. The solid line corre- 
sponds to q = 0.2, the long-dashed line to the standard model with q = 1, and the 
short-dashed line to q = 5. 

5 Discussion 

The time-dependent hydrodynamic simulations presented in this paper clearly 
show that, within the limitations of the idealized model, tidally- induced spiral 
shocks can lead to significant radial transport in accretion disks in binary 
systems. Our results agree qualitatively with previous numerical work; namely, 
we find that tidal forces excite a two-armed spiral shock wave that remains 
stationary in the rotating frame of the binary system. However, our use of high 
spatial resolution and a high-order numerical algorithm that conserves angular 
momentum about the accreting star allows us to quantitatively measure the 
transport effects of these spiral waves. 

Perhaps the most interesting result is that a is large. In all of our simulations, 
the effective a was at least 0.1 in the outer regions of the accretion disk 
where the tidal shock waves are excited. However, as predicted by analytic 
arguments, a decays with decreasing radius. In cold disks such as those found 
in CV systems, this drop is so precipitous that we could not measure an a 
below a radius of only 1/3 that of the outer edge of the disk. 

These numerical results agree in many respects with the analytic work of Lar- 
son (1989) and Spruit (1987). The shape of the spiral wave agrees reasonably 
well with the simple prediction that the trailing angle should be roughly the 
inverse of the orbital Mach number (Godon 1997), and the Mach number of 
the spiral shocks matches the values found in the self-similar models of Spruit 
(1987). Even the value of the effective a in terms of the Mach number of the 
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spiral shocks, M s , agrees with the analytic estimate of Larson (1989). 

However, the value of the effective a in our simulations is much larger than 
the analytic predictions of Larson (1989) and Spruit (1987) when written in 
terms of the orbital Mach number of the disk. In addition, the dependence of 
a on sound speed and disk radius does not agree with their analytic theory, 
which predicts an effective a proportional to (c;?r) 3 / 4 . This discrepancy thus 
appears to lie in the derivation of the Mach number of the spiral shock as a 
function of the sound speed in the disk. 

We note, however, that the analytic work applies to the propagation of spiral 
waves through a disk, and does not give reference to the source of the waves. 
This work shows that tidal forces drive a strong two-armed spiral wave, which 
in turn produces an effective a « 0.1 in the outer regions of the disk. As 
this wave propagates inward, there is a competition between growth and dis- 
sipation, as well as the effects of an increasing Mach number. In our models 
dissipation wins out and the wave weakens (as measured by radial velocity, 
density compression, or effective a) as it propagates to smaller radii. Perhaps 
it would eventually approach the analytic solution as it continues to weaken, 
but the corresponding a of such a wave is much smaller than we would be 
able to measure. 

A further result of these simulations is the relative independence of spiral 
waves on the mass ratio in the binary system. The strength of the two-armed 
spiral shocks at their origin near the outer edge of the disk was fairly constant 
in all of our models. This has the consequence that, despite all else, one can 
be confident that the effective a in the outer regions of an accretion disk in a 
binary system is (at least) ~ 0.1. 

Finally, these models allow us to investigate the effects of tidal truncation. 
We found a strong dependence of the tidal radius on the Mach number in 
the disk, such that colder disks extended out to larger radii. The tidal radius 
in our coldest model (with q = 1) agrees well with the analysis of Paczynski 
(1977), who found the largest non- intersecting periodic orbit in a pressure-less 
disk. Our results suggest that pressure effects in a warm disk can significantly 
decrease the tidal radius (30% smaller for a Mach number of ~ 10 at the outer 
edge of the disk). 

There are three important caveats concerning our idealized model of an accre- 
tion disk. The first is our assumption of an isothermal equation of state. In- 
ternal heating in a real accretion disk will lead to a sound speed that increases 
with decreasing radii. Spiral waves in such a disk will not be as hampered by 
the rapid increase in Mach number that we found in our simulations, so we 
might expect the transport due to such waves to be significant deeper into the 
disk. A real equation of state (with 7 > 1) is also expected to produce stronger 
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transport effects (Spruit 1987). Our second major assumption is the limita- 
tion to two dimensions. In general, one might expect the inclusion of the third 
dimension to result in weaker spiral waves in the equatorial plane, since the 
vertical direction provides for the possibility of additional wave modes, vertical 
expansion that may diminish the effects of shock compression, or refraction of 
waves out of the equatorial plane. 

The third caveat is the missing tidal stream in our model. We found the 
most important effects of the spiral shocks to occur at the outer edge of the 
accretion disk, yet for a steady-state accretion disk there must be some source 
of incoming gas. In a close binary system this is most likely a tidal stream 
originating at the inner Lagrange point of the binary potential. This stream 
will presumably strike the outer edge of the disk and drive its own spiral wave 
into the disk (Shu 1976). Will this impact wave transport angular momentum? 
Will it interfere with the tidally-driven spiral waves? 

While there is still a great deal of work to be done, it is clear from this 
and previous work that tidally-driven spiral waves can effectively drive mass 
transport through the outer regions of an accretion disk in a close binary 
system. 
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